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CHAPTER  1 


MATHEMATICAL  OVERVIEW 


1.1  Introduction 

Program  AGAUS  treats  electromagnetic  scattering  and  absorption  of  spherical 
particles  by  Mle  theory.  The  rigorous  formulation  of  the  problem  of  the 
Interaction  of  a  sphere  and  a  plane  wave  can  be  found  In  a  number  of  books  on 
advanced  electromagnetic  theory1  2  3  and  will  not  be  repeated  here.  In  the 
Interest  of  clarity  and  uniformity  of  notation,  however,  a  brief  summary  of 
the  needed  results  of  the  theory  is  Included.  The  primary  objective  of  this 
section  is  to  provide  descriptions  of  ways  in  which  AGAUS  handles  integrations 
over  distributions  of  particle  sizes  and  mixtures  of  different  particle  size- 
distributions. 

In  analytic  notation,  a  "size-distribution"  (or,  a  "model")  is  described  by  a 
"distribution  function"  f(r)  with  the  property  that  f(r)dr  is  the  relative 
number  of  particles  (within  a  unit  volume  of  space)  having  radii  between  r  and 
r  +  dr.  (In  a  strict  mathematical  sense,  f(r)  is  a  "probability  density 
function,"  but  the  present  usage  is  common.)  The  total  number  of  particles 
per  unit  volume  (N)  is  found  by  integrating  f(r)dr  between  minimum  and  maximum 
values  of  r  appropriate  to  a  given  physical  situation.  Within  AGAUS,  f(r)  has 
the  units  "particles  per  cubic  centimeter  per  micrometer  radius."* 

1.2  Notations  for  Individual  Particles 

Common  Mfe  theory  assumes  that  the  Interacting  sphere  is  homogeneous  and 
characterized  by  a  radius  and  a  complex  index  of  refraction  n  =  m  -  ik.  The 
overall  interaction  is  usually  characterized  by  quantities  called  "efficiency 
factors"  represented  by  the  letter  Q.  Four  such  quantities  are  used  in  AGAUS 
and  are  given  subscripts  e,  a,  s,  and  r  representing  extinction,  absorption, 
scattering,  and  back  scattering,  respectively.  When  multiplied  by  the 
geometric  area  of  the  sphere,  the  efficiency  factors  give  cross  sections  for 
the  above  processes.  Cross  sections  are  often  represented  by  the  letter  C, 
with  an  appropriate  subscript.  The  Q's  and  C's  are  functions  of  the  complex 
index  of  refraction  (itself  a  function  of  wavelength)  and  of  the  radius  of  the 
sphere  and  the  wavelength  of  incident  radiation.  Except  for  the  implicit 


Wan  de  Hulst,  1957,  Light  Scattering  by  Spherical  Particles,  John  Wiley  and 
Sons,  Inc.,  New  York 

2W.  K.  H.  Panofsky  and  M.  Phillips,  1962,  Classical  Electromagnetic  Theory, 
second  edition,  Addison-Wesley  Publishing  Co.,  Cambridge,  MA,  pp  233 

3M.  Born  and  E.  Wolf,  1959,  Principles  of  Optics,  Pergamon  Press,  Oxford, 
England,  pp  630 

♦With  the  exception  of  type  0,  within  the  actual  computer  code,  f(r)  is 
internally  normalized  in  such  a  way  that  N  =  1,  and  the  absolute  number  of 
particles  per  unit  volume  is  determined  from  data  that  must  be  supplied  by  a 
user.  For  this  reason,  f(r)  is  called  a  "relative"  number  density  in  chapter 
2. 
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dependence  of  m  and  k  on  the  wavelength,  the  particle  radius  and  wavelength 
affect  the  efficiency  factors  only  In  ratio  form,  giving  rise  to  a 
dimensionless  “size  parameter"  a— the  circumference  of  the  sphere  divided  by 
the  wavelength. 

The  detailed  angular  distribution  of  scattered  radiation  is  described  by  a 
"phase  function"  p(n),  which  gives  the  scattered  energy  per  unit  solid  angle 
per  unit  incident  radiation.  The  phase  function,  like  the  Q's,  is  a  function 
of  radius,  the  wavelength,  and  the  index  of  refraction.  Since  problems 
treated  by  A6AUS  are  azimuthal ly  symmetric,  p  is  often  written  as  p(e)  even 
though  it  Is  properly  a  function  of  both  e  and  $. 

1.3  Notations  for  a  Single  Size-Distribution  Model 

In  dealing  with  an  aerosol  model  consisting  of  many  particles  of  identical 
composition  but  varying  radii,  the  efficiency  factors  (and  the  phase  function) 
must  be  weighted  by  the  distribution  function  f(r)  to  secure  an  overall  cross 
section.  Thus,  for  example,  the  total  scattering  cross  section  for  all 
particles  in  a  unit  volume  of  space  is  given  by 


CSU) 


*r2Qs(r,m,k,x)  f(r)dr 


(1.1) 


with  similar  formulae  for  extinction,  etc. 

In  the  same  vein,  one  can  define  an  average  cross  section  per  particle  (or  a 
cross  section  per  average  particle)  via 


Zs  =  CSU)/N  , 


where  N  is  the  number  of  particles  per  unit  volume, 
function  p(e)  is  found  from 


(1.2) 


The  average  phase 


.  max 

p(e)  =  jj-  /  p(Q)  f(r)dr  . 


(1.3) 


1.4  Mixtures  of  Size-Distribution  Models 

To  indicate  how  AGAUS  combines  mixtures  of  two  or  more  aerosol  models,  a 
subscript  1  will  be  added  to  f(r)  and  to  the  Q's  and  p.  Thus  f^(r)  will 

represent  the  distribution  function  for  the  1th  component  of  a  mixture,  etc. 
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(1.4) 


Then 


max 

N.  =  J  f.  (r)dr  , 
! 

rmi  n 


C  (X)  =  /irr2Q  (r,m,k)f  .(r)dr 
s.  s1  1 


(1.5) 


Fe  =  Cc  (X )/N. ,  etc. 
S1  S1  1 


(1.6) 


The  total  cross  sections  for  a  mixture  are  then  found  from  formulae  such  as 


-  I 

i 


"l*s. 


(1.7) 


The  combined  phase  function  is  calculated  from 


p(a) 


{  “i  Pi(9)  rSl 


*  "l  cs 


i 


(1.8) 


so  that 


/  /  p  (ii )  dsi  =  • 


(1.9) 


Note  that  the  FINAL  phase  functions  printed  by  A6AUS  are  normalized  such  that 
the  above  Integral  is  unity. 

The  detailed  relationships  between  the  Q's,  p(e),  and  the  parameters  m,  k,  and 
a  are  summarized  in  section  1.8. 

1.5  Extinction,  Scattering,  and  Radar  Coefficients 

Readers  may  note  that  the  "cross  sections"  discussed  above  do  not  really  have 
the  usual  dimension  of  area,  but  rather  the  dimensions  of  area  per  unit 
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volume,  or  of  inverse  length.  Thus,  for  example,  Ce  represents  the  extinction 

cross  section  for  all  those  particles  located  in  a  unit  volume  of  space. 
Multiplication  of  Ce  by  the  volume  of  an  optical  path  is  required  to  get  the 

true  extinction  cross  section  for  the  path.  If  one  assumes  a  unit  cross 

section  for  such  a  path,  then  Ce  can  be  regarded  as  an  extinction 

"coefficient"  Kg.  AGAUS  calculates  and  prints  several  such  "coefficients" 

having  the  final  unit  "per  kilometer."  As  an  example,  consider  a  homogeneous 

path  of  L  kilometers.  The  optical  depth  of  the  path  is  given  by 


*  »  Ke  L  , 


(1.10) 


and  the  radiation  intensity  I ( L )  would  be  found  from 


I(L)  =  1(0)  exp(-x) 


(1.11) 


1.6  Treatment  of  Hygroscopic  Particles 

Hygroscopic  particles  absorb  water  from  an  ambient  atmosphere  and  may  undergo 
growth  to  sizes  that  differ  from  their  initial  values.  The  easiest  type  of 
growth  to  formulate  is  that  In  which  the  absorbed  liquid  water  and  the  initial 
material  form  a  solution  droplet.  Even  this  simple  type  of  situation  is  quite 
complicated,  but  certain  useful  approximations  do  exist  and  are  used  in 
AGAUS.  The  best  known  survey  of  the  growth  of  hygroscopic  particles  and 
attendant  changes  in  optical  properties  Is  that  of  G.  Hanel,^  whose  reading  is 
highly  recommended  for  anyone  wishing  to  utilize  these  features  of  AGAUS. 

The  starting  point  for  size  adjustments  is  the  equation4 


1  + 


lV(aw) 


(1.12) 


4G.  Hanel,  1976,  "The  Propertle*  of  Atmos  herlc  Aerosol  Particles  as  Functions 
of  the  Relative  Humidity  at  ermodyr  .c  Equilibrium  with  the  Surrounding 
Moist  Air,"  Advances  in  Geophysics  Vo  19,  Academic  Press,  New  York 


which  is  essentially  a  combined  form  of  equations  (2.31)  and  (2.22)  of 
reference  4.  The  quantities  appearing  here  are  ps,  the  mass  density  of  the 

original  aerosol  material  (the  solute);  pw,  the  mass  density  of  pure  liquid 

water;  »,  Hanel's  "growth  factor";  and  aw,  the  "water  activity"  of  the  system. 

The  water  activity  aw  is  related  to  the  equilibrium  particle  radius  r,  the 

fractional  humidity  f  (percent  per  100),  and  the  temperature  T  through  the 
equation 


(1.13) 


,  2ov 

where  c  =  , 

in  which  a  is  the  molecular  weight,  v  Is  the  surface  tension,  R  is  the 
universal  gas  constant,  T  is  the  kelvin  temperature. 

Equations  (1.12)  and  (1.13)  are  coupled  by  r  and  aw  and  are  solved 
simultaneously  by  means  of  an  iterative  procedure. 

From  equation  (1.13) 


(1.14) 


Substitution  of  equation  (1.14)  into  (1.12)  yields 


(1.15) 


Equation  (1.15)  is  then  used  as  follows;  an  Initial  estimate  of  r  Is  taken  to 
be  the  dry  radius  r  and  is  used  to  calculate  the  right-hand  side  of  the 
equation.  This  procedure  provides  a  new  value  of  r  via  the  equality.  The  new 
value  is  then  used  to  recalculate  the  right-hand  side  of  equation  (1.15)  to 
get  a  second  estimate  of  r.  This  procedure  Is  continued  until  a  new  estimate 
differs  from  the  previous  one  by  no  more  than  0.01  percent. 


« 

4G.  Hanel,  1976,  "The  Properties  of  Atmospheric  Aerosol  Particles  as  Functions 
of  the  Relative  Humidity  at  Thermodynamic  Equilibrium  with  the  Surrounding 
Moist  Air,"  Advances  In  Geophysics,  Vol  19,  Academic  Press,  New  York 
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After  exit  from  the  Iteration  procedure,  equation  (1.13)  Is  then  used  to 
calculate  a  value  for  aw.  The  net  droplet  density  is  then  found  from 


where 


“s'1 


uaw 
1  -  a. 


-)/A 


(1.16) 


(1.17) 


Equations  (1.16)  and  (1.17)  are  essentially  equation  (2.46)  of  reference  4. 
The  complex  index  of  refraction  n  *  m  -  Ik  Is  then  obtained  from 


and 


m  =  mw  +  (m$  -  mJ/A  , 


(1.18) 


k  =  kw  +  (k$  -  kw)/A  . 


(1.19) 


Since  the  growth  factor  u  Is  generally  a  function  of  aw,  the  procedure 

outlined  above  is  not  truly  complete.  A  user  must  still  supply  a  value  for  u 
In  some  fashion.  As  a  matter  of  fact,  y  may  even  be  an  indirect  function  of 
radius,  but  Hanel  states  that  It  Is  reasonable  to  use  an  average  value  u  and 
to  Ignore  that  complexity. 

Perhaps  the  easiest  way  to  enter  a  value  of  y  into  A6AUS  Is  via  the  input 
parameter  EMUA  (17—  assumed  to  represent  an  averaged  value  as  noted  above)  and 
to  Ignore  the  dependence  of  y  on  r.  EMUA  must  still  be  regarded  as  a  function 
of  aw  (or  the  fractional  humidity),  however.  The  only  reliable  way  to 

determine  the  appropriate  values  of  7  for  a  given  material  Is  by  experimental 

measurements.  Reference  4  provides  some  guidance  In  the  matter  as  well  as 
some  examples  of  y  -  f  dependencies  for  several  aerosol  models. 


4G.  Hanel,  1976,  "The  Properties  of  Atmospheric  Aerosol  Particles  as  Functions 
of  the  Relative  Humidity  at  Thermodynamic  Equilibrium  with  the  Surrounding 
Moist  Air,"  Advances  In  Geophysics,  Vol  19,  Academic  Press,  New  York 


to 


1.7  Validation  of  AfiAUS 


The  continued  fraction  subroutine  (MIEGX),  which  does  the  actual  Mie 
calculations  in  AGAUS,  was  evaluated  by  Burlbaw  and  Miller5  through 
comparisons  with  J.  V.  Dave's  routine  DAMIE.6  In  all  cases  for  which  the 
latter  routine  was  considered  to  be  reliable  (size  parameters  to  400  for 
nonabsorbing  spheres,  and  smaller  size  parameters  for  absorbing  spheres),  the 
agreements  were  to  five  digits  or  better.  MIEGX  was  also  tested  against  a 
code  MIE2,  developed  by  Radiation  Research  Associates,  and  again  showed 
equally  good  agreement.  Unlike  most  other  known  codes,  MIEGX  does  not  fail  at 
the  machine  level  (that  is,  does  not  produce  underflows  or  overflows  arising 
from  finite  precision  "roundoffs,"  etc)  even  at  size  parameters  of  1000  or 
greater  (see  also  section  2.2.3). 

The  halving  method  of  integration  employed  in  AGAUS  has  been  validated  by 
comparing  the  results  of  complete  runs  for  several  different  aerosol  models 
with  results  published  by  others,  including  Shettle  and  Fenn7  and 
Dei rmendjian.B  Some  of  these  evaluations  have  been  reproduced  in  reference 
5.  Any  users  who  wish  to  make  similar  evaluations  should  remember  principles 
of  the  halving  method*  and  should  be  sure  to  make  the  convergence  criterion 
(delta)  small  enough  to  insure  that  the  minimum  interval  in  radius  is 
comparable  to  that  used  by  others. 

1.8  Sumary  of  Mie  Theory  Used  In  AfiAUS 

Mie  theory  predicts  the  scattering  by  and  the  absorption  in  an  isolated, 
discrete,  homogeneous,  isotropic  sphere  of  radius  r  with  a  known  complex 
refractive  index  n  =  m  -  ik  relative  to  the  surrounding  medium  and  illuminated 
by  monochromatic  radiant  energy  with  wavelength  x  in  the  surrounding  medium. 
The  theory  is  given  in  detail  in  standard  texts  and  need  not  be  repeated 
here.  Instead,  those  elements  of  theory  needed  for  an  understanding  of  the 
numerical  algorithms  used  are  included.  The  information  in  this  section  is 
partially  based  on  material  from  Burlbaw  and  Miller5  and  Shirkey  et  al.9 


5E.  Burlbaw  and  A.  Miller,  1981,  Modification  of  Single  Scattering  Model 
AGAUSX,  ASL-CR-0780-1 ,  US  Army  Atmospheric  Sciences  Laboratory,  White  Sands 
Missile  Range,  NM 

6J.  V.  Dave,  1972,  Development  of  Programs  for  Computing  Characteristics  of 
Ultraviolet  Radiation^  Technical  Report  for  Contract  NAS5-21860  ( tJASA) ,  IBM 
Corporation 

7E.  P.  Shettle  and  R.  W.  Fenn,  1979,  Models  for  the  Aerosols  of  the  Lower 
Atmosphere  and  Effects  of  Humidity  Variations  on  their  Optical  Properties, 
AFGL-TR-79-0214,  Air  Force  Geophysics  Laboratory,  Hanscomb  Field,  MA 

BD.  Deirmendjian,  1969,  Electromagnetic  Scattering  on  Spherical 

Polydispersions,  American  Elsevier  Publishing  Company,  Inc,  New  York 

*See  section  1.10. 

9R.  C.  Shirkey  et  al,  1930,  Single  Scattering  Code  AGAUSX:  Theory, 

Applications,  and  Listing,  ASL- 1  k-uuo*,  (JS  Wmf  Atmospheric  sciences 
Laboratory,  White  Sands  Missile  Range,  NM 


All  scattering  properties  of  spheres  are  computed  from  m  and  k,  and  through 
the  use  of  the  induced  electric  and  magnetic  multipole  moments  of  the  sphere 
at  and  b£,  respectively.  The  moments  are  given  by* 


fj(na)  T£(a)  -  n  ?£(na)  f£(a) 

"  tj(na)  C4(oJ  -  n  T£lna)  ^(a)  ’  (1.20) 

and 

n  f^not)  *£(a)  -  TA(na )  ^(o) 
bt  =  n  Y'(na)  Cjjla)  -  f£lna)  S£(a)  *  (1.21) 

The  prime  denotes  differentiation  with  respect  to  the  argument.  The  v£(z)  and 

5£(z)  functions  are  Ricatti -Bessel  functions  of  the  first  and  third  kinds, 

respectively,  and  are  related  to  spherical  Bessel  functions  j  (z)  and  n  (z)  by 

£>  £ 


T£(Z)  =  Z  j£(z)  , 


(1.22) 


and 


C£(z)  =  zj£(z)  -  1  zn£(z)  =  H'£(z)  +  1Xjl(z)  , 


(1.23) 


where 


jt(z) 


>/2 


J1+1/2*Z^  • 


(1.24) 


*a  Is  the  Mie  size  parameter  and  is  equal  to  2nr/x,  where  r  is  the  radius  of  a 
sphere. 
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and 


(1.29) 


Although  the  cross  sections  account  for  the  energy  removed  from  the  forward 
beam,  they  do  not  give  any  information  about  where  the  scattered  photons  go. 
This  Information  is  contained  in  scattering  amplitudes  and  intensity  factors 
that  relate  the  flux  density  scattered  through  an  angle  a  relative  to  the 
incident  flux  density.  There  are  two  amplitudes,  S T ( c )  and  S2 ( e ) ,  and 

intensity  factors  i x ( e)  and  i 2 ( e ) ,  which  correspond  to  1  ight,  respectively, 

polarized  perpendicular  and  parallel  to  the  plane  of  scattering  defined  by  the 
direction  of  incidence  and  the  direction  of  scattering. 

The  intensity  factors  are  related  to  the  scattering  amplitudes  by 


i  i  ( 0 )  =  |  S  i  ( 9 )  | 2  , 
i2(9)  •  |S2(0)|2  . 


The  amplitudes  come  from  the  multi  pole  moments  through 


s>(e)  ■ i£4Vv«(9>  ♦  v«(9»  - 


*=i 


(1.32) 


and 


s*(9)  =  XctKVv.19**  vt(e)1  • 


(1.33) 


and  angular  factors  *  (0)  and  t  (0)  defined  in  terms  of  associated  Legendre 
functions: 


nA(0)  =  P*  (cos  0)/sin  0 


(1.34) 


dPi (cos  0) 


Alternative  expressions  frequently  used  are 


* t (0)  = 


dP. (cos  0) 


(1.36) 


dw  (e) 

*»<«>  =  cos  6  •  ’i(e>  -  s1"2  8  •  WcoTi 


(1.37) 


where 


P  (cos  0)  =  - - (cos2  0  -  l)1  . 

1  2lt\  d(cos*  0) 


These  functions  satisfy  the  following  recurrence  relations: 


(1.38) 


,t(e)  •  cos  .  .  (.)  - 


t-hrv2(e)  • 


(1.39) 


t£{0)  =  COS  0  Ut(e)  "  ■  ^sin20  ff4-l^e^  +  n)l-2(9) 


(1.40) 


The  scattering  cross  section  measures  the  ability  of  a  particle  to  scatter 
light;  consequently,  Csca  is  expected  to  be  obtained  from  an  integral  over  the 

scattering  intensity  factors.  Equation  (1.27)  follows  from 


Csca  "  77  [i  +  i2(0))  d(cos  0)  *  (1.41) 
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Although  the  Intensity  factors  themselves  may  be  used  In  scattering 
calculations,  they  are  primarily  suited  for  computing  flux  densities,  and 
frequently  scattered  light  can  be  more  conveniently  measured  and  computed  In 
terms  of  radiances.  Radiances  do  not  have  an  inverse  square  distance 
dependence;  therefore,  the  distance  from  the  scatterer  to  the  detector  need 
not  be  known  If  the  detector  field  of  view  is  small  and  Is  filled  by  the 
scattering  cloud.  The  phase  function  p(e)  gives  a  radiance  I  scattered  Into 
the  9  direction  in  terms  of  the  radiance  I0  Incident  on  the  particle. 

The  phase  function  is  dimensionless  and  is  defined  in  AGAUS  as 


p(  e  >  = 


Lii(e)  +  i2(e)]  . 


(1.42) 


The  integral  of  the  final  normalized  phase  function  printed  by  AGAUS  is  unity. 

For  the  special  case  of  e  =  180°,  backscatter,  the  efficiency  is  expressed  by 
the  radar  cross  section  o.  The  radar  cross  section  may  be  defined  as  4*  times 
the  back  scattered  power  per  steradlan  divided  by  the  incident  power  per  unit 
area  or 


o  *  4irr2I(r,180°)/Io  . 


(1.43) 


This  expression  can  be  reduced  by  the  relations 


I(r,e) 


I0Hi(e)  +  1 2  ( © )  > 

2k2  r2 


(1.44) 


i  i (180° )  =  12(180°)  =  (S^IBO0)!2  . 


(1.45) 


o  -  ZL  JS1(180°)|2  , 


(1.46) 
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and  when  divided  by  the  geometrical  cross  section,  G  =  nr2, 


4|Si  (180°)p 

0  -  0  a  _________ 

wradar  E  02 


(1.47) 


where  a  =  2*r/A  .  Using 


-*t(180°)  =  t£(180°)  =  (-1)1  •  \  l  (t  +  1)  , 


(1.48) 


one  obtains 


-S^ISO0)  =  S2(180°) 


oo 

l  («  +7)  (-l)V  -  b  )  . 

i*l  L  l  l 


(1.49) 


1.9  Numerical  Techniques  for  Mle  Theory  In  AGAUS 

Most  of  the  numerical  calculations  performed  by  AGAUS  are  those  associated 
with  calculating  the  coefficients  a£  and  b^.  These  calculations  are  performed 

by  a  subroutine  (MIEGX)  that  is  based  on  a  "continued- fraction"  method  for 
finding  the  Ricatti -Bessel  functions.  The  method  employed  is  that  of  W.  J. 
Lentz.10  This  method  was  chosen  because  It  provides  very  good  accuracy  using 
only  single-precision  arithmetic  and  because  of  the  ease  with  which  it  handles 
large  size-parameters  and  imaginary  parts  of  the  complex  index  of 
refraction.  Since  the  method  is  not  widely  employed,  it  is  reviewed  in  some 
detail  here. 

Subroutine  MIEGX  used  In  AGAUS  Is  a  modified  version  of  a  routine  (DAMIE) 
originally  developed  by  J.  V.  Dave.6  The  major  modification,  made  by  Lentz,10 
was  the  use  of  a  continued-fraction  method  rather  than  ordinary  backward 
recursion  techniques  for  evaluating  some  of  the  quantities  appearing  in 
equations  (1.20)  and  (1.21).  The  continued-fraction  method  will  be  discussed 
In  more  detail  below.  It  is  noted,  at  this  point,  however,  that  the  angular 
functions  n  and  t  are  evaluated  in  MIEGX  using  forward  recursion  and 


10W.  J.  Lentz,  1976,  "Generating  Bessel  Functions  In  Mle  Scattering 
Calculations  using  Continued  Fractions,"  Appl  Opt,  15(3) :668 

6J.  V.  Dave,  1972,  Development  of  Programs  for  Computing  Characteristics  of 
Ultraviolet  Radiation  technical  Report  for  contract  NASb-2i»bU  (NASA),  18m 
Corporation 


equations  (1.39)  and  (1.40)  and  starting  values  r  (e)  =  0,  wjfe)  =  1, 
t © ( 0 )  *  0,  and  x^e)  =  cos  e  .  o 

For  computational  purposes  of  MIEGX,  equations  (1.20)  and  (1.21)  may  be 
rewritten  In  the  following  forms:6 


•»t(no) 

Re[c£(a)]  -  Re[?£_^(a)] 

.  n  a. 

*4 

A  (na)  * 

4  +  4 

n  aj 

54(a)  •  C4-l(a) 

(1.50) 


b 


4 


nA.(na)  +  — 

4  a 

■ 

nA.(na)  +  — 
k  a 


ReC^fa)]  -  ReU^ta)] 
C£(a)  -  C£-1(a) 


(1.51) 


where 


(1.52) 


and  z  -  no  is  complex. 

Several  methods  have  been  used  for  evaluating  the  quantities  5£  and  A^.  The 
principal  difference  between  MIEGX  and  Dave's  Mle  routines  is  the  way  in  which 
the  ratios  A^  are  found.  Lentz's  continued-fraction  method  evaluates  the 
quantities  A£.  A  brief  summary  of  the  latter  method  . Is  given  below,  but 
Interested  readers  should  refer  to  Lentz's  article  for  detailed  Information. 

To  simplify  notation,  the  common  form  for  a  continued  fraction: 


6J.  V.  Dave,  1972,  Development  of  Programs  for  Computing  Characteristics  of 
Ultraviolet  Radiation,  Technical  Report  for  Contract  NAS5-21860  (NASA),  IBM 
Corporation 
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will  be  written  as 


f  -  91©-g-©-H-©*3-  +  •••  • 

1  *2  3  4 


(1.55) 


With  the  first  of  those  notations,  Lentz  showed  that 


J£-l(z)  24 rs 

ozr'“e 


-© 


i 


-© 


i 


2  ~  _i  ~  _i 

(4  +  l)z  2(4  +  2 ) z  2(4  +  3 ) z 


T  +  **•  ‘  (1.56) 


To  clarify  the  procedure  used  in  M1EGX,  it  is  most  convenient  to  use  a  second 
compacted  notation  for  the  continued  fraction  (1.54): 


f  =  [9  •  9  »  9  »  •  •  • 3  • 

12  3 

Then,  defining  the  "kth  convergent"  to  be 

fk  "  ^i»  V  9 3 . » 

Lentz  showed  that 

Cg  !»••• »9  3C9l ,9  ] 

1  *_a  1  *  1 

f|c  * - . 

[9^].. •  [g^_^ » ••  •  »9  Kg^ , . .. ,g ^3 


(1.57) 


(1.58) 


(1.59) 
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The  order  of  the  terms  appearing  in  equation  (1.59)  has  been  reversed  to 
emphasize  that  the  evaluation  begins  with  g}  rather  than  with  an  unknown  value 
9k* 

The  ratio  J  .(z)/J,(z)  can  now  be  expressed  in  the  form  of  equation  (1.59)  if 
one  takes 


gi  =  (-l)i+12U  +  i  +  1)  z'1 


(1.60) 


Subroutine  MIEGX  utilizes  the  foregoing  procedures  for  calculation  of  the 


ratios  J, 


_l/J ^  needed  to  evaluate  the  Mie  coefficients  a^  and  b^. 


practice,  MIEGX  initially  determines  the  value  of  and  some  value  of  t,  say 

i  =  N,  such  that  N  is  the  largest  value  of  t  whose  use  is  expected  in  a 

particular  run.  The  required  ratio  Is  evaluated  by  calculating  the 

relevant  value  of  fk  [equation  (1.59)]  for  a  sequence  of  k  values  up  to  the  k 
for  which 


1]  <  e  , 


(1.61) 


where  c  is  a  tolerance  value  corresponding  to  the  precision  available  within  a 

particular  computer.  e  is  given  the  value  of  10~6  for  32-bit  precision 
calculations. 

Once  An  has  been  found,*  the  values  of  A^  for  i  <  N  are  found  by  backward 
recursion  and  are  stored  in  an  array  for  recall  as  needed  in  evaluating  the 
Mie  summations  of 


Re[S  (0)],  Im[S  (e)],  Re[S  (e)]  and  Im[S  (e)] 
112  2 


The  sum  is  terminated  when 


*If  N  exceeds  the  dimension  NOIM  of  the  complex  array  A(  ),  MIEGX,  then 
a2xNDIM  (°r  a3xNDIM*  etc*)  *s  calculated  and  backward  recursion  is  used  to 

find  the  needed  A.*, 
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Kl2  *  |btl2 


(1.62) 


and  when  the  fractional  change  In  the  radar  efficiency  is  also  less  than  e, 
that  is,  when 


^rad  t  "  *Vad  t-1 


<  e  . 


(1.63) 


This  is  more  stringent  than  the  first  test  alone  and  is  a  test  on  the  phase 
functions  as  well. 

For  size  parameters  smaller  than  0.01,  the  cumbersome  Mie  apparatus  is 
bypassed  and  Van  de  Hulst's1  analytic  approximations  are  used  for  finding  the 
efficiency  factors  and  phase  function. 

MIEGX  returns  the  following  quantities  as  required  by  AGAUS:  Qext»  ^sca* 

Qrad,  and  P(J).  The  P(J)  as  returned  by  MIEGX  are  average  intensities  and 

must  be  further  normalized  to  become  the  actual  phase  functions  [equation 
(1.42)].  Some  limitations  associated  with  numerical  Implementation  of  the 
above  procedures  are  discussed  in  section  2.2.3. 

1.10  Other  Numerical  Procedures  Used  In  AGAUS 

AGAUS  handles  numerical  Integrations  of  the  type  indicated  in  equations  (1.3) 
and  (1.5)  using  the  "trapezoidal  rule"  and  a  process  called  "halving."  The 
latter  process  is  described  next. 

A  general  Idea  of  how  program  AGAUS  performs  integrations  over  sizes  can  be 
obtained  by  considering  a  numerical  method  for  determining  the  area  under  a 
curve  g(r)  versus  r  such  as  that  sketched  in  figure  1.  The  objective  is  to 
evaluate  numerically  an  Integral 


max 

G  =  /  g(r)dr 


(1.64) 


^an  de  Hulst,  1957,  Light  Scattering 
Sons,  Inc.,  New  York  ' 


>her1cal  Particles,  John  Wiley  and 
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to  some  desired  degree  of  accuracy  using  the  smallest  number  of  values  of  r 
for  the  numerical  calculations.  The  procedure  adopted  in  AGAUS  is  as  follows 
(reference  to  figure  1  may  be  helpful): 

a.  An  initial  estimate  Gj  of  the  value  of  the  integral  is  made  with  three 
values  of  r  labeled  by  Roman  numeral  1  and  the  "trapezoidal  rule." 

b.  A  second  estimate  G2  is  then  made  with  increments  Ar,  which  are  half 

as  large  as  those  used  in  getting  Gj .  In  getting  G2,  the  two  additional  r- 
values  labeled  II  are  utilized. 

c.  The  values  of  G2  and  G  are  compared  to  each  other  by  calculating  a 
quantity 


|  G2  -  Gx  | 

6  *  ~1^T“ 


(1.65) 


and  comparing  it  to  a  preset  quantity  a. 

If  6  <  a,  then  it  is  assumed  that  G2  is  a  "sufficiently"  accurate 

representation  of  G,  and  the  computations  are  terminated.  If  on  the  other 
hand,  6  >  a,  one  proceeds. 


g(r) 


Figure  1.  Sketch  of  how  values  of  particle  radii  are  selected  in  the  halving 
procedure  of  AGAUS. 
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d.  A  third  estimate  G3  of  G  is  then  made  by  again  cutting  the  increment 

Ar  to  half  its  previous  value.  This  results  in  the  addition  of  confutations 
at  the  three  new  r-values  labeled  III  in  figure  1.  A  new  value  of  6  is  then 
calculated  from 


|G3  -  G2 


(1.66) 


and  6  is  again  compared  to  a.  If  6  is  still  greater  than  a,  the  spacing 
between  r-values  is  cut  in  half  once  more,  and  the  "estinati on-comparison" 
process  is  repeated  until  either  6  <  a  or  some  maximum  number  of  r-values  has 
been  reached. 

In  AGAUS,  the  quantity  used  in  the  testing  process  is  the  total  volume  of  the 
aerosol  particles,  namely, 


V  *  /g-  *r3  f(r)dr 


(1.67) 


Volume  was  chosen  as  the  quantity  to  be  tested  because  its  r3  dependence  might 
make  it  converge  more  slowly  than  the  extinction  or  scattering  cross  sections 
while  maintaining  a  reasonable  convergence  for  phase  functions.  The  code  may 
be  easily  modified  to  test  on  any  desired  quantity;  however,  the  maximum 
number  of  radii  allowed  by  AGAUS  depends  on  the  number  of  size  groups  defined 
in- ‘a  run  and  may  be  changed  by  revision  of  the  source  code  (see  section  2.3.2 
for  additional  discussion  of  this  matter). 


CHAPTER  2 
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2.1  Introduction 

AGAUS  is  a  Mie  scattering  code  for  calculation  of  the  single-scattering 
properties  of  polydisperse  spherical  particles.  The  primary  quantities 
calculated  oy  AGAUS  are  the  extinction,  scattering,  absorption,  and 
backscatter  coefficients  plus  the  angular  intensity  distribution  produced  by 
various  forms  of  particle  size-distributions.  The  code  contains  seven 
different  size-distribution  models  and  also  permits  a  user  to  provide  an 
arbitrary  size  model  as  a  set  of  discrete  numerical  input  data.  There  is,  in 
addition,  a  special  subcode  that  may  be  used  to  tabulate  the  Mie  efficiency 
factors  (extinction,  scattering,  absorption,  and  backscatter)  as  a  function  of 
particle  radius.  An  important  feature  of  this  version  of  AGAUS  is  the  ability 
to  treat  and  combine  the  scattering  properties  of  a  mixture  of  up  to  five 
different  size-distributions  in  a  single  run  and  to  repeat  the  process  for  an 
arbitrarily  large  number  of  wavelengths.  A  method  for  handling  the  effect  of 
relative  humidity  on  the  equilibrium  sizes  of  hygroscopic  particles  is  also 
i ncluded. 

2.2  Stnaary  of  Output  Quantities 

This  section  presents  a  summary  of  the  quantities  produced  and  printed  in  a 
run  of  AGAUS.  More  detailed  descriptions  may  be  found  in  subsequent  sections. 

2.2.1  Definitions 

a.  "Size-distribution  model"  refers  to  a  particular  functional 
relationship  between  the  relative  number  of  particles  per  unit  radius  and  the 
values  of  the  radius. 

b.  "Di stribution-type"  refers  to  a  number  that  identifies  the  various 
"size-distribution  models"  that  are  recognized  by  AGAUS  as  valid. 

c.  "Size-group"  or  "size-bin"  refers  to  particular  ranges  of  values  of 
particle  radii  for  which  certain  optical  properties  are  calculated  and  printed 
by  the  program.  To  clarify  what  is  meant,  consider  a  situation  in  which  one 
wants  to  examine  particles  with  radii  between  O.Oum  and  10um,  but  would  like 
to  see  the  properties  of  particles  whose  radii  lie  in  the  ranges  0.0pm  to 
2.0pm,  2.0pm  to  5.0pm,  and  5.0pm  to  10.0pm.  This  would  be  handled  by  dividing 
the  overall  range  (0.0  to  10.0)  into  three  "size-groups,"  with  "group  1" 
belonging  to  radii  between  0.0pm  and  2.0pm,  etc. 

2.2.2  Major  Outputs 

With  those  definitions  the  major  outputs  can  now  be  described: 

a.  Particles  per  Gran.  The  number  of  particles  in  a  1-g  sample  of 
material  if  all  the  particles  in  the  sample  had  radii  belonging  to  a  size- 
group. 
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b.  Haber  Fraction.  The  ratio  of  the  number  of  particles  that  would 

belong  to  a  size-group  to  the  sun  of  the  numbers  of  particles  in  all  size- 

groups  being  used  In  a  run.  As  an  example,  supDOse  that  three  size-groups  are 

being  used  and  that  a  1-g  sample  contains  10  particles  with  sizes  belonging  to 
group  1,  100  particles  in  group  2,  and  5  particles  in  group  3.  Then  the 
number  fraction  for  group  1  would  be  10/(10  +  100  ♦  5).  For  group  2  it  would 
be  100/(10  +  100  +  5),  and  so  forth. 

c.  Mass  Fraction.  The  fraction  of  the  mass  contained  by  particles  in  ALL 

size-groups  that  would  belong  to  particles  whose  radii  place  them  in  a 
selected  size-group.  Continuing  the  above  example,  suppose  that  the  10 
particles  in  group  1  has  a  mass  of  0.0005  g,  those  in  group  2  a  mass  of  0.80 
g,  and  those  in  group  3  a  mass  of  0.1995  g.  Then  the  mass  fraction  for  group 

1  would  be  0.0005;  for  group  2,  0.8;  and  for  group  3,  0.1995. 

d.  Cumulative  Mass  Fraction.  The  sum  of  the  mass  fractions  of  all  size- 
groups  up  to  and  including  the  group  for  which  this  quantity  Is  generated. 
The  cumulative  mass  fractions  for  the  examples  above  would  then  be  0.0005, 
0.8005,  and  1.0000. 

e.  Group  Cross  Sections  In  Square  Centimeters  per  Gram.  Extinction, 
absorption,  scattering  cross  sections  that  would  exist  if  1  g  of  material 
existed  among  particles  whose  radii  fell  inside  a  size-group. 

f.  Group  Cross  Sections  In  Square  Centimeters  per  Particle.  The 

quantities  of  the  preceding  form  divided  by  the  number  of  particles  per  gram 
(i tern  1  above) . 

g.  Weighted  Averages.  (Extinction,  absorption,  etc.)  cross  sections  as 
integrated  against  the  relative  number  of  particles  per  unit  radius  described 
by  the  size-distribution  model  in  use.  These  results  have  the  units  of  square 
centimeters  per  particle. 

h.  Extinction,  Scattering,  and  Back scattering  Coefficients.  The  weighted 
averages  (Item  7,  above)  multiplied  by  the  total  number  of  particles  per  cubic 
centimeter  associated  with  a  particular  size-model  after  summation  over  all 
size-distribution  models  being  treated  In  a  particular  run  of  AGAUS.  The  unit 
of  these  quantities  Is  "per  kilometer,"  so  that  the  "optical  depth"  of  a 
transmission  path  is  the  product  of  these  "coefficients"  and  the  length  of  the 
path  in  kilometers. 

i.  Phase  Functions.  The  relative  intensity  per  unit  solid  angle  of 
scattered  radiation  as  a  function  of  angle  as  averaged  over  all  pertinent 
size-groups  and  size-distribution  models.  The  phase  functions  are  normalized 
so  that  their  integrals  over  all  solid  angles  are  unity. 

j.  Attenuation  Coefficient.  This  quantity  gives  the  total  extinction 
coefficient  with  units  of  square  meters  per  milligram  of  aerosol  material  to 
be  found  In  an  optical  path.  The  attenuation  coefficient  may  appear  with  two 
“suffixes"  referring  to  "dry"  and  "wet"  aerosols.  In  that  context,  "dry" 
means  the  number  of  milligrams  of  aerosol  present  before  any  adjustment  of 
particle  radii  caused  by  accretion  of  liquid  water  arising  from  the 
hygroscopic  properties  of  the  material.  "Wet"  means  the  number  of  milligrams 
of  the  aerosol  that  would  exist  after  taking  hygroscopic  properties  into 
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account.  In  other  words,  the  “wet"  attenuation  is  that  which  would  be 
expected  after  the  dispersal  of  1  mg  of  a  hygroscopic  material  into  an 
environment  whose  relative  humidity  will  lead  to  the  accretion  of  liquid  water 
and  an  Increase  in  the  sizes  and  masses  of  the  original  aerosol  particles 
(that  Is,  when  Hanel's  growth  formulae  are  invoked). 

In  addition  to  the  above  quantities,  A6AUS  also  provides  a  printed 
recapitulation  of  most  of  the  input  data  and  several  other  derived  results 
that  may  be  of  interest  to  at  least  some  users.  If  warranted,  certain 
diagnostic  and  warning  messages  are  also  printed. 

2.3  Input  data  Requirements 

Input  data  for  A6AUS  are  formulated  as  a  set  of  standard  80-column  punched 
cards  or  card-images.  At  least  four  types  of  data  cards  are  required  for  a 
run.  If  a  user  selects  certain  program  options,  as  many  as  three  additional 
types  of  data  cards  may  be  needed.  Those  optional  data  cards  are  identified 

below  as  card  types  3A,  3B,  and  3C.  Again  depending  on  options  selected,  some 

types  of  cards  may  need  to  be  repeated  to  provide  additional  data  such  as  a 
new  wavelength  value. 

2.3.1  Card  Type  1 

This  card  is  required  (although  it  may  be  blank)  and  may  contain  up  to  80 
alphanumeric  characters  Identifying  the  run.  The  format  is  40A2. 

2.3.2  Card  Type  2 

This  card  contains  a  set  of  II  program  control  parameters  entered  as  integers 
with  the  format  (1115).  The  symbolic  names  of  these  parameters  and  their 
meanings  are  described  next  and  nwst  appear  on  the  data  card  from  left  to 
right  in  the  order  in  which  they  are  given  below. 

NWAVE  NWAVE  is  nominally  the  number  of  wavelengths  to  be  treated  within  a 

single  run  of  the  program,  but  may  also  be  used  to  define  the  number 
of  times  one  wishes  to  loop  over  any  variable  found  on  cards  of  type 
5  (below).  NWAVE  must  be  1  or  greater.  A  minimum  of  NWAVE  cards  of 
type  5  is  necessary. 

NIDSTP  NIDSTP  is  the  number  of  different  size-distribution  types  to  be 
combined  at  each  wavelength.  The  maximum  value  of  NIDSTP  Is  5. 

IW  IW  is  a  parameter  that  may  be  used  to  Indicate  that  the  particles  of 

all  size-distribution  models  are  to  be  treated  as  though  they  were 
composed  of  liquid  water.  To  select  that  option,  IW  should  be  set  to 
zero.  If  IW  =  0  any  values  of  the  complex  refractive  index  and  the 
mass  density  (specific  gravity)  found  on  card  type  5  will  be  Ignored 
and  replaced  by  those  of  liquid  water  at  the  temperature  (°C)  that  is 
taken  from  card  type  5.  To  avoid  these  features,  set  IW  =  1.  If  IW 
=  0,  the  wavelength  on  data  card(s)  type  5  mist  be  between  0.2pm  and 
200.0pm. 
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IDELT  IOELT  is  used  to  signal  the  presence  or  absence  of  card  3A  from  the 
input  data  deck.  If  card  type  3A  is  not  to  be  present,  set  IDELT  to 
zero.  Use  IDELT  =  1  if  card  type  3A  will  be  present.  If  IDELT  =  0, 
default  values  of  the  data  that  may  appear  on  card  3A  are  defined  by 
the  program. 

NANG  Parameter  NANG  specifies  the  number  of  angles  at  which  the  phase 
function  is  to  be  calculated.  The  value  of  NANG  oust  be  between  1 
and  65.  If  phase  functions  are  of  no  particular  interest,  set  NANG  = 
1  or  0.  The  values  of  the  angles  used  in  the  calculation  will  depend 
on  the  parameter  IANG  (next).  The  value  of  NANG  has  a  significant 
effect  on  running  time,  so  one  should  use  as  small  a  value  as  is 
consistent  with  other  needs. 

IANG  IANG  is  used  to  choose  from  among  three  ways  of  selecting  the  angles 
at  which  phase  functions  are  calculated.  Note,  however,  that  the 
phase  functions  are  printed  only  if  parameter  NANG  is  3  or  greater. 

IANG  =  0  will  override  the  input  value  of  NANG  (setting  it  to  65)  and 
will  result  in  the  use  of  65  equally  spaced  angles  between  0  and  180 
degrees. 

IANG  ■  1  will  recognize  the  value  of  NANG  as  entered  and  result  in 
the  use  of  NANG  equally  spaced  angles  between  0  and  180  degrees. 

IANG  =  2  informs  the  program  that  cards  of  type  3B  carrying  values 
for  the  angles  to  be  used  will  be  found  in  the  data  deck.  This 
option  allows  a  user  to  specify  whatever  set  of  values  for  the  angles 
may  be  desired. 

NBINS  NBINS  is,  in  general,  the  number  of  size-groups  into  which  the  run  is 
to  be  subdivided.  (NBINS  =  0  is  an  exception  discussed  later.)  The 
maximum  permissible  value  of  NBINS  is  50.  If  NBINS  is  not  zero,  a 
user  must  supply  cards  of  type  3C  in  the  input  data  deck.  Those 

cards  carry  the  user's  definitions  of  the  radii  to  be  associated  with 
each  size-group.  An  input  value  of  zero  for  NBINS  will  provide  the 
use  of  just  one  size-group,  with  the  minimum  and  maximum  radii 
determined  either  from  data  entered  on  cards  of  type  4  or  by  the 
program  itself.  If  NBINS  *  0,  there  must  be  no  data  cards  of  type  3C 
in  the  input  data  deck.  Note  also  that  NBINS  must  be  set  to  zero  if 
use  of  size-distribution  type  number  zero  is  selected  on  a  card  of 
type  4  below.  (See  also  section  2.3.2.) 

IEO  Parameter  IEO  controls  the  generation  of  an  output  data  disc  or  tape 
file  containing  several  run  parameters  (below)  and  the  values  of  the 
phase  functions.  If  no  such  file  is  wanted,  be  sure  to  set  IEO  = 
0.  There  are  two  other  legal  choices  for  the  value  of  IEO: 

IEO  *  1  overrides  the  users  values  of  NANG  and  IANG  and  results 
in  the  use  of  65  predetermined  angles  for  phase  function 
calculations.  This  option  will  result  in  a  file  compatible  with 
EOSAEL  usage. 
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IEO  *  2  also  overrides  the  value  IANG  on  card  2  and  sets  IANG  = 
2.  The  latter  means  that  the  user  mist  supply  a  set  of  data 
cards  of  type  38  carrying  NANG  values  for  the  angles  to  be  used 
in  the  phase  function  calculations. 

In  both  cases  the  output  file  is  written  on  logical  unit  number 
NEOU  (next  parameter).  The  file  will  contain  the  following 
quantities: 

1.  The  angles,  with  format  (11(F6.2,1X) ) . 

2.  The  number  of  angles,  a  phase  function  identifier  (set 

to  zero  by  AGAUS,  which  Implies  user  supplied  data  for 

subsequent  use  as  EOSAEL  input),  the  wavelength 
(micrometers),  the  albedo  for  single  scattering,  and  the 
extinction  and  scattering  coefficients  (per  kilometer). 
These  items  are  written  with  the  format  (2(12, IX), 
F5.2,1X,F8.6,1X,2(E12.6,1X)). 

3.  The  phase  functions,  written  with  a  format  of 
(6(E12.6,1X)),  at  the  angles  noted  immediately  above.  The 
phase  functions  are  normalized  so  that  their  integrals  over 
all  solid  angles  are  unity. 

NEOU  NEOU  is  the  logical  unit  number  to  be  used  to  identify  the  file 

controlled  by  IEO.  If  In  doubt,  use  zero  for  both  IEO  and  NEOU. 

NUNIT  If  NUNIT  is  nonzero,  the  following  quantities  are  written  on  logical 
unit  number  NUNIT  for  each  radius  at  which  Mie  calculations  are  made 
during  a  run:  dry  radius,  radius  after  any  Increase  related  to 

relative  humidity,  the  relative  particle  number  at  the  dry  radius, 
and  the  Mie  efficiency  factors  for  extinction,  scattering  and 

back scatter. 

MQRTE  Normally  zero,  but  if  MQRTE  =  12345,  the  quantities  listed  under 

NUNIT  will  be  directed  to  whatever  logical  unit  has  been  assigned  to 
the  listing  device  for  the  run.  (That  assignment  is  logical  unit 
number  6  unless  altered  by  changing  the  value  of  the  parameter  IOUT 
in  the  COMMON  BLOCK  subroutine  of  AGAUS). 


These  cards  are  more  or  less  optional  but  their  presence  or  absence  in  a  data 
deck  must  be  properly  correlated  with  the  values  of  the  parameters  IOELT, 
IANG,  IEO,  and  NBINS,  which  appear  on  card  2. 


Card  3A  contains  two  variables:  DELTA  and  ATOP  with  the  format  (2E10.4). 


DELTA  DELTA  is  the  value  of  the  "convergence  criterion"  used  to  decide  when 
to  terminate  the  Integration  over  particle  sizes  within  each  size- 
group.  The  running  time  of  AGAUS  Is  Intimately  connected  to  the 
value  of  DELTA  used  In  a  run.  The  default  value  of  DELTA  (used  If 
IDELT  *  0  and  card  3A  is  therefore  absent)  Is  0.001.  A  value  of  0.01 
will  probably  meet  most  needs  If  great  accuracy  Is  not  needed.  A 
value  of  0.1  is  often  adequate  for  sample  or  test  runs. 


ATOP 


This  parameter  sets  an  upper  limit  to  the  Hie  size-parameter  for 
which  complete  Mie  calculations  are  done.  Its  default  value  is 
400.  The  function  of  ATOP  is  the  following:  whenever  a  particle 
radius  has  a  value  that  corresponds  to  a  si ze-parameter  that  is 
larger  than  ATOP,  the  actual  value  of  the  size-parameter  is  ignored 
and  a  size- parameter  equal  to  ATOP  is  used  instead.  The  first  time 
that  happens,  the  results  of  the  Mie  calculation  are  saved  and 
substituted  when  any  subsequent  particle's  size-parameter  exceeds 
ATOP.  This  procedure  can  save  computation  time  and  usually  does  not 
introduce  a  significant  amount  of  error  in  the  values  of  the 
efficiency  factors  (and  cross  sections)  because  they  approach 
constant  values  at  large  size-parameter  values.  The  effect  of  the 
approximation  is  more  severe  on  the  phase  functions,  however,  but  are 
not  likely  to  be  devastating  since  the  relative  mmber  of  very  large 
particles  is  ordinarily  small.  If  ATOP  is  made  larger  than  400, 
inherent  numerical  behavior  of  the  Mie  subroutine  may  lead  to 
unreliable  values  of  the  backscatteri ng  efficiency  factor  and 
possibly  of  the  phase  functions  associated  with  large  size- 
parameters. 

The  actual  value  of  ATOP  at  which  the  above-mentioned  instabilities 
begin  to  appear  is  linked  to  the  size  of  a  complex  array  A(  )  in 
subroutine  MIEGX  rather  than  to  a  more  intrinsic  property  such  as  use 
of  single  precision  arithmetic.  Those  instabilities  first  seem  to 
appear  when  the  size  parameter  being  treated  reaches  a  value  that  is 
approximately  equal  to  the  dimension  of  the  array  A.  In  the  standard 
version  of  MIEGX,  A  is  defined  to  contain  400  complex  elements  or  800 
single  precision  floating  point  numbers.  The  choice  of  400  rather 
than  200  or  1000  was  made  as  a  compromise  between  program  core 
requirements  and  anticipated  Mie  size  parameters  found  in  normal 
usage.  If  sufficient  memory  is  available,  the  dimension  of  A(  )  may 
be  increased  substantially,  thereby  moving  the  size  parameter  at 
which  instabilities  might  occur  to  values  larger  than  400.  Test  runs 
have  been  made  with  the  dimension  of  A  as  large  as  3000  without  the 
appearance  of  obvious  instabilities  of  the  type  under  discussion. 
Users  who  contemplate  increasing  the  dimension  of  A(  )  must  note  that 
one  other  parameter  (NDIM)  would  also  require  a  change.  This  is 
discussed  explicitly  in  paragraph  e  of  section  2.4.3. 

Card(s)  3B  carry  the  angles  to  be  used  in  the  phase  function  calculations  if  a 
user  chooses  to  supply  them.  These  cards  carry  up  to  16  values  per  card  with 
the  format  (16F5.1).  The  symbolic  quantity  read  from  these  cards  is  an  array 
labeled  ANGL  in  the  source  code.  The  number  of  entries  should  be  equal  to  the 
value  of  NANG  placed  on  card  2.  There  must  be  no  extra  cards:  for  NANG 
between  i  and  16,  use  one  card;  for  NANG  between  17  and  32,  use  two  cards, 
etc. 

Card(s)  3C  are  used  to  supply  the  upper  radius  limit  for  the  size-groups, 
which  are  read  under  the  format  (8E10.4).  The  symbolic  quantities  on  this 
type  of  card  are  labeled  RR  (a  51-element  array)  in  the  FORTRAN  source  code. 
The  values  of  the  RR's  must  be  In  ascending  order,  and  the  first  entry  should 
be  the  upper  limit  in  radius  for  the  first  size-group.  (The  lower  limit  on 
radius  for  the  first  size  group  is  normally  zero  if  these  cards  are  used  at 
all,  that  is,  if  NBINS  is  nonzero).  There  must  be  at  least  NBINS  entries,  and 


no  extra  cards  are  permitted.  There  must  of  course  be  enough  cards  to  make 
the  number  of  data  values  equal  to  N8INS.  Remember  that  NBINS  nust  not  be 
greater  than  50. 

2.3.4  Card  Type  4 

Type  4  cards  carry  the  data  that  select  the  size-distribution  model  (or 
models)  to  be  assumed  for  a  run  and  certain  parameters  associated  with  the 
different  models.  Up  to  five  cards  of  this  type  are  permitted  for  most  runs 
(exception  is  type  0;  see  below).  The  number  of  cards  of  this  type  that  are 
required  is  specified  by  the  value  of  NIDSTP  on  card  2. 

All  cards  (except  some  cards  for  model  IDSTP  =  0;  see  below)  have  the  general 
form: 


IDSTP,  Ql,  Q2,  Q3,  Q4,  Q4,  Q6 

The  format  is  (I3,7X,6E10.4)  and  in  all  cases  IDSTP  is  an  integer  that 
identifies  a  "distribution-type."  It  must  be  placed  in  column  3.  The  Q's 
have  different  interpretations  for  different  distribution-types  and  will  be 
given  symbolic  names  in  the  ensuing  lines. 

2.3.4. 1  Type  0.  Arbitrary  User-Supplied  Discrete  Model.  This  distribution- 
type  allows  a  user  to  enter  from  l  to  !>13  pairs  of  data  (one  pair  per  card) 
for  the  value  of  a  particle  radius  and  an  absolute  number  of  particles  having 
that  radius.  This  model  is  selected  by  using  a  type  4  data  card  containing  a 
zero  in  column  3  and  the  number  of  data  pairs  as  a  floating  point  number  in 
columns  11-20.  In  other  words,  Ql  (above)  on  the  first  type  4  card  should  be 
the  number  of  data  pairs  to  be  expected  by  the  program.  Note  that  this 
particular  distribution  type  requires  an  absolute  number  density. 

A  user  must  then  insert  NRADI  data  cards  each  carrying  a  value  of  radius  (in 
micrometers)  and  an  absolute  number  density  (in  particles  per  cubic  centimeter 

per  micrometer),  where  NRADI  Is  the  number  of  radius  values  to  be  used.  The 

format  is  (2E10.4). 

Users  should  probably  not  try  to  mix  a  type  0  model  with  some  other  type  under 
the  NIDSTP  greater  than  1  option.  Such  mixtures  have  not  been  tested. 

2.3. 4. 2  Type  1.  Lognormal  Distribution 

IDSTP,  RL0,  RHI,  RBAR,  SIGMA,  DENS 

Ql  is  RL0  -  The  minimum  radius  at  which  the  model  is  to  be  cutoff. 

Q2  is  RHI  -  The  maximum  radius  at  which  the  model  is  to  be  cutoff. 

Q3  is  RBAR  -  The  radius  at  which  the  relative  number  density  peaks. 

Q4  Is  SIGMA  -  The  standard  deviation  to  be  used  (Not  LN(SIGMA)). 

Q5  is  DENS  -  The  number  of  particles  per  cubic  centimeter. 


31 


This  distribution  has  the  mathematical  form: 


f(r)  =  ± 


exp 


1 

7 


|£n(r/r)| 
e  n<  c ) 


(2.1) 


where  f(r)  dr  is  the  relative  number  of  particles  with  radii  between 
r  and  r  +  dr,  and  C  is  a  normalization  constant. 

RLO  and  RHI  are,  in  a  sense,  optional:  If  RLO  =  RHI,  then  the  code 
will  Itself  assign  cutoff  values  (that  is,  will  redefine  both)  that 
insure  that  most  of  the  particle  radii  will  be  included  in  the 
calculation.  To  select  that  option  just  make  the  values  of  Q1  and  Q2 
equal  to  zero.  Note,  however,  that  they  nust  not  be  omitted:  RBAR, 
SIGMA,  and  DENS  must  still  be  in  positions  Q3,  Q4,  and  Q5.  This 
option  relative  to  RLO  and  RHI  is  not  valid  for  the  other 
distribution  types  (below). 

DENS  is  also  optional  but  only  in  a  sense.  If  a  user  does  not  want 
to  specify  the  number  density,  entering  a  value  of  0.0  will  tell  the 
program  to  calculate  the  number  density  from  the  mass  density  (RHOA, 
card  5),  the  mass  concentration  (CONC,  card  5)  and  the  (computed) 
average  volume  per  particle.  If  both  DENS  and  CONC  are  zero,  the 
program  assigns  a  default  value  DENS  -  1.0.  This  behavior  of  DENS  is 
also  applicable  to  the  other  distribution  types  (listed  below)  that 
also  use  DENS  as  an  input  parameter. 

2. 3.4. 3  Type  2.  Double  Exponential  Model 

IDSTP,  RLO,  RHI,  Q,  A,  B,  DENS 

RLO,  RHI,  and  DENS  have  the  same  meaning  as  listed  under  type  1. 

The  meanings  of  A,  B,  and  Q  are  best  described  by  the  mathematical 
form  for  the  distribution: 


f (r)  =  Qe‘A/r  +  (1  -  Q)e'B/r 


(2.2) 


(In  the  FORTRAN  source  code  Q  has  the  symbolic  name  CUE.) 
2.3. 4. 4  Type  3.  Power  Law  Model 

IDSTP,  RLO,  RHI,  Q,  A,  DENS 
The  form  of  the  distribution  is 


f (r)  =  Qr 


RHO,  RHI,  and  DENS  are  again  (In  fact,  are  always)  defined  as  in  type 
1  above. 


2. 3. 4. 7  Type  6.  Modified 


Fog  Model 


IOSTP,  RLO,  RHI,  RC,  ALF,  GAM,  ELWC 


The  form  of  this  model  Is  the  same  as  in  type  5.  This  model  always 
assumes  that  the  aerosol  material  will  be  liquid  water  droplets  and 
finds  the  value  of  OENS  from  the  value  of  ELWC.  ELWC  Itself  is  the 
"liquid  water  content"  In  units  of  grams  per  cubic  meter.  Use  of 
this  model  relieves  a  user  from  the  need  to  specify  the  optical 
constants  (EMA  and  CAYA,  card  5)  and  the  mass  density  (RHOA,  card  5) 
and  an  Internal  table  lookup  procedure  provides  those  data  for 
water.  The  tables  of  Hale  and  Querry  (Applied  Optics,  12:555,  ff) 
are  used  in  this  procedure. 


2. 3.4.8  Type  7.  Special  Table  Generation  Mode.  This  is  not  actually  a 
distribution  model  at  all.  Out  is  an  operating  model  used  for  calculating  and 
printing  taoles  of  Mie  efficiency  factors  as  a  function  of  radius.  The  input 
data  are  as  follows: 


IDSTP,  RIO,  RHI,  OELR 

RLO  is  the  minimum  particle  radius  to  be  used. 
RHI  is  the  maximum  particle  radius  to  be  used. 
l)ELR  is  the  increment  in  radius  to  be  used. 


This  "distribution-model,"  like  all  the  others,  will  expect  the  user 
to  provide  the  wavelength  and  the  optical  constants  at  each 
wavelength  on  data  card  type  5. 

The  quantities  that  are  printed  are  the  adjusted  (per  Hanel)  radius, 
the  dry  radius,  the  corresponding  Mie  size-parameter,  the  efficiency 
factors  for  extinction,  scattering,  absorption  and  backscatter,  and 
the  adjusted  (again  per  Hanel)  index  of  refraction.  The  program 
formats  50  sets  of  results  per  page. 

This  "mode"  may  be  used  to  generate  several  sets  of  tables  in  a 
single  run,  but  doing  so  requires  a  nonstandard  way  of  arranging  the 
input  data  cards  of  types  4  and  5.  NWAVE  on  card  type  2  should  be  1; 
NIDSTP  (card  1)  should  be  equal  to  the  number  of  "sets"  for  the  run 
(but  no  greater  than  5);  NIDSTP  pairs  of  data  cards  of  types  4  and  5 
must  be  "interleaved"  rather  than  collected  together  as  distinct 
groups  of  NIDSTP  type  4  cards  followed  by  NWAVE  type  5  cards. 

2. 3.4.9  Type  8.  Marshall -Palmer  Rain  Model 

IDSTP,  RAIN 


RAIN  is  the  rainrate  in  millimeters  per  hour. 


The  form  for  f(r)  is 


f(r)  =  Ae'^r 


(2.7) 


with  A  =  1.6E-5  and  B  =  8.2E-4  *  (RAIN  **  (-0.21)).  The  value  of 
DENS  is  found  from  DENS  =  A/B,  which  is  theoretically  valid  only  if 
RLO  =  0  and  RHI  -  infinity.  In  practice  the  distribution  is  cutoff 
at  RLO  =  0.0001  and  RHI  =  2500.0  (micrometers),  but  this  should 
induce  no  appreciable  error  in  DENS  since  few  rain  drops  will  exceed 
2.5  mm  In  radius. 

A  few  words  of  caution  may  be  helpful  in  connection  with  use  of  this 
model.  The  maximum  Mie  size- parameter  used  is  roughly  13,000  divided 
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by  the  wavelength.  Since  the  Mie  subroutine  used  in  AGAUS  becomes 
unreliable  at  size-parameters  of  the  order  of  400  to  1000  (depending 
on  which  quantities  are  of  most  interest),  the  model  is  only  useful 
at  fairly  large  wavelengths  and  is  intended  for  use  at  wavelengths  of 
the  order  of  1  mm  or  so.  Because  the  data  tables  for  optical 
constants  for  water  do  not  extend  beyond  200um,  use  of  the  Marshal  1- 
Palmer  model  is  generally  inconsistent  with  setting  IW  =  0  on  data 
card  2,  and  users  must  provide  the  appropriate  optical  constants  on 
cards  of  type  5  (below). 

2.3.4.10  Type  9.  Hybrid  Lognormal  and  Power-Law  Model 

IDSTP,  RBAR,  SIGMA,  RMIN,  RMAX,  PWR 

This  model  consists  of  the  lognormal  model  (type  1)  for  radii  between 
zero  and  RMIN  and  a  power-law  (type  3)  model  for  radii  between  RMIN 
and  RMAX.  The  number  density  functions  f(r)  are  made  equal  at  the 
radius  defined  by  the  value  of  RMIN. 

RBAR  and  SIGMA  have  the  same  definitions  as  in  type  1. 

PWR  is  the  same  quantity  as  A  In  type  3. 

If  either  RBAR  or  SIGMA  Is  zero,  this  model  becomes  the  same  as  type 
3,  with  RLO  =  RMIN  and  RHI  =  RMAX.  If  RMIN  and  RMAX  are  equal,  then 
the  model  is  the  same  as  type  1. 

One  restriction  Imposed  by  computer  exponent  limits  is  known  to  exist 
for  the  hybrid  model,  namely  that  RMIN,  RBAR,  and  SIOIA  must  be  such 
that  the  quantity 


I £n(  RMIN/RBAR) i 

|-TnT5raU)  1 


13  . 


2.3.5  Card(s)  Type  5 

Cards  of  this  type  carry  data  describing  the  wavelength  to  be  used  and  the 
general  physical  properties  of  the  aerosol  materials  whose  relative  number 
distributions  are  specified  by  cards  of  type  4.  At  least  one  of  the  type  5 
cards  is  required  for  each  run  of  AGAUS.  More  generally  the  number  of  type  5 
cards  needed  is  given  by  the  product  of  NWAVE  and  NIDSTP.  The  format  of  type 
5  cards  Is  (8E10.4).  The  symbolic  names  for  the  data  are: 

WAVE,  EMA,  CAYA,  RHOA,  CONC,  RELHUM,  TEMP,  EMUA 

WAVE  -  the  wavelength  In  micrometers 

EMA  -  the  real  part  of  the  index  of  refraction  of  the  aerosol 


CAYA  -  the  absolute  value  of  the  Imaginary  part  of  the  index  of 
refraction 
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RHOA  -  the  mass  density  (specific  gravity)  of  the  aerosol  in  grams 
per  cubic  centimeter 

CONC  -  the  mass  concentration  (grams  of  aerosol  material  per  cubic 
centimeter  of  space) 

RELHUM  -  the  relative  humidity  in  percent 

TEMP  -  the  environmental  temperature  in  Celsius  degrees 

EMUA  -  Hanel ' s  growth  factor  "mu-bar" 

One  card  of  this  type  must  be  supplied  for  each  value  of  NWAVE  and  for  each 
value  of  IDSTP  at  each  wavelength,  regardless  of  the  value  of  IW  on  card  2. 

2.4  Additional  Program  Information 

2.4.1  Choosing  Values  for  DELTA 

As  was  stated  under  the  description  for  data  card  type  3,  DELTA  controls  the 
exit  from  the  numerical  integration  routine  in  which  the  efficiency  factors 
are  converted  to  cross  sections  and  weighted  by  the  relative  number  density, 
etc.  The  integration  is  done  using  the  trapezoidal  rule  in  a  series  of  steps 
called  "halvings."  Each  size-group  is  treated  separately.  On  the  first  pass 
across  a  size-group,  computations  are  done  at  the  endpoints  and  midpoint  of 
the  range  of  radii  belonging  to  the  group.  In  addition  to  the  cross  sections, 
weighted  sums  of  the  phase  functions  and  the  total  volume  of  the  particles  are 
computed.  The  program  then  calculates  new  estimates  of  all  the  above 
quantities  by  performing  the  computations  at  radii  that  lie  halfway  between 
the  radii  used  in  the  first  pass  and  combines  the  newest  results  to  form  new 
estimates  of  the  integrals.  After  the  first  two  passes,  five  values  of  radius 
will  have  been  used.  The  proqram  then  proceeds  through  a  series  of  loops  in 
which  the  spacing  of  the  radii  is  made  half  those  used  in  the  most  recent 
pass.  At  the  end  of  each  pass  (after  the  second),  the  most  recent  value  for 
the  integrated  particle  volume  is  compared  to  the  one  that  existed  at  the  end 
of  the  previous  pass.  Exit  from  this  procedure  occurs  as  soon  as  the 
difference  between  the  "new"  and  "old"  values  for  volume  divided  by  the  "new" 
value  is  less  than  DELTA.  If  this  criterion  is  not  satisfied  after  some 
preset  number  of  "halvings"  have  occurred,  the  program  exits  anyway  and  prints 
a  message  to  the  effect  that  convergence  was  not  achieved  in  treating  that 
particular  interval.  The  number  of  radii  used  before  automatic  exit  is 
determined  by  an  internal  program  variable  NHALV  (accessible  in  the  source 
code  but  not  as  a  normal  input  parameter),  and  its  value  is  dependent  on  the 
number  of  size-groups  being  treated  in  a  given  run. 

The  maximum  number  of  radii  used  for  each  size-group  is  as  follows: 


If  NBINS  is  0,  1,  or  2,  513  radii  may  be  used  before  exit. 

If  NBINS  is  3  or  4,  only  257  radii  per  group  are  allowed. 

If  NBINS  is  5  or  larger,  129  radii  are  allowed  per  size-group. 


These  assignments  can  be  changed  In  the  source  code’s  main  program  if  more  or 
fewer  radii  are  desired,  except  for  the  type  zero  model.  Models  other  than 
type  zero  calculate  the  new  radius  values  as  needed,  out  type  zero  sets  a 
maximum  of  513  radii  and  places  those  values  in  an  array  for  recall.  This 
procedure  is  used  because  the  user  supplied  type  zero  data  are  unlikely  to  be 
spaced  at  the  correct  values  of  radius  required  by  the  halving  process  and  an 
interpolation  is  made  from  the  user  supplied  values  of  the  radii  to  those 
needed  for  the  halving  procedure.  Users  who  have  ample  memory  available  could 
increase  the  value  of  513  to  1025,  2049,  etc.,  by  redefining  the  dimensions  of 
the  arrays  R{  ),  F{  ),  and  FFF(  )  in  the  source  code  and  recoding  a  portion  of 
the  interpolation  procedure.  If  this  process  is  contemplated,  pay  particular 
attention  to  the  parameter  NHALV:  the  maximum  number  of  radii  must  be  equal 
to  1  +  (2  raised  to  the  power  NHALV) . 

The  running  time  of  the  program  is  mainly  set  by  the  number  of  times  the  Mie 
subroutine  is  called  and  the  values  of  the  size-parameters  in  use;  therefore, 
the  use  of  a  very  small  value  of  DELTA  will  force  the  program  to  utilize  the 
maximum  allowed  number  of  radii  for  each  size-group.  The  number  of  radii  used 
might  well  be  far  greater  than  that  needed  for  acceptable  precision  in  the 
results  and  hence  make  the  run  require  an  excessive  amount  of  time.  The 
situation  is  most  critical  when  the  number  of  size-groups  is  fairly  large.  In 
many  cases  the  overall  contributions  of  the  particles  in  some  size-groups  to 
the  total  extinction  may  be  negligible  because  of  the  small  relative  number  of 
particles  in  those  groups,  and  use  of  a  very  small  value  of  DELTA  can  result 
in  a  considerable  amount  of  computation  time  for  large  particles  that  often  do 
not  really  affect  the  results  very  much.  The  program  does  not  (with  one  type 
of  exception)  make  any  judgments:  the  same  value  of  DELTA  is  used  for  all 
size-groups.  The  exception  is  as  follows:  in  runs  with  NBINS  greater  than 
zero,  the  program  finds  a  quantity  called  VNIAX  for  each  size-distribution 
model  in  use.  This  parameter  consists  of  the  total  particle  volume  (weighted 
by  occurrence  frequency)  of  particles  for  whatever  size-group  makes  that 
quantity  a  maximum  for  the  whole  range  of  groups  in  use.  VMAX  is  used  in  the 
integration  loops  discussed  above  in  such  a  way  that  an  automatic  exit  from 

computations  for  a  size  group  will  occur  if  the  volume  estimated  for  the  first 

33  choices  of  radius  is  less  than  VMAX*l.E-6.  This  procedure  can  easily  be 
disabled,  if  a  user  so  desires,  by  removing  a  single  line  from  the  source  code 
for  subroutine  AGXP2.  The  line  in  question  is  one  that  reads 

IF(N.GT.5.AND. VOLHHT.LT. VMAX) GOTO  ... 

Suggested  values  for  DELTA  are  in  the  range  0.01  to  0.001.  If  after  choosing 
values  of  DELTA  in  this  range,  the  program  outputs  look  "strange"  or  appear  to 
be  questionable,  then  a  rerun  with  a  smaller  value  of  DELTA  may  be 

advisable.  Users  who  have  not  previously  used  AGAUS  should,  perhaps,  do  a 
little  experimenting  with  a  few  different  values  of  DELTA  to  see  how  its  value 
affects  their  results. 

For  distribution  type  zero  (arbitrary,  user-supplied  model),  the  igested 
value  of  DELTA  is  0.001  or  smaller  because  the  interpolations  used  . n  that 
model  can  present  problems  if  DELTA  is  not  fairly  small. 


2.4.2  Selecting  Values  for  NBINS  and  the 
Size-Group  Radius  Llelts 


In  many  Instances,  users  will  not  really  care  to  subdivide  the  range  of  radii 
to  be  used  in  a  run  into  more  than  one  size-group.  In  such  cases,  it  will 
probably  be  sufficient  to  set  the  parameter  NBINS  =  0,  and  to  control  the 
maximum  and  minimum  cutoff  radii  via  the  parameters  RLO  and  RHI  on  data  cards 
of  type  4  (or,  for  IDSTP  =  9,  via  RMAX) .  Some  situations  have  arisen, 
however,  in  which  choosing  NBINS  greater  than  zero  is  advisable  even  though 
the  breakdown  of  outputs  by  size-groups  is  not  particularly  wanted.  An 
example  of  such  a  situation  is  one  in  which  the  NIDSTP  >  0  option  is  being 
invoked  to  treat  a  bimodal  (NIDSTP  =  2)  or  multimodal  model  in  which  the 
distributions  that  are  being  combined  peak  at  different  values  of  radius.  In 
that  kind  of  situation,  it  is  advisable  to  break  the  computation  range  into 
intervals  that  are  roughly  centered  on  the  radii  at  which  the  various  modes 
have  their  peaks.  Doing  that  will  help  the  user  in  being  assured  that  enough 
points  are  treated  for  each  mode  to  make  the  convergence  checks  meaningful. 

To  clarify  this  state  of  affairs,  consider  an  aerosol  model  being  described  by 
the  lognormal  distribution  IDSTP  =  1,  with  the  first  mode  being  peaked  at  r  = 
0.05um  and  a  standard  deviation  of  about  1.5.  Suppose  the  other  mode  peaks  at 
5.0pm  and  has  a  standard  deviation  of  2.0.  To  assure  coverage  of  the  major 
part  of  the  second  mode,  one  would  probably  choose  RLO  to  be  something  of  the 
order  of  0.001  and  RHI  to  be  of  the  order  of  100.  Now,  a  run  with  NBINS  =  0 
permits  a  maximum  of  513  radii  to  be  used  in  the  Mie  calculations  and 
Integrations  over  sizes.  The  Increment  between  adjacent  radii  (if  all  allowed 
radii  are  used  in  the  halving  process)  is  (RHI-RL0)/513,  or  about  0.2pm.  This 
increment  Is  far  too  large  to  give  adequate  attention  to  the  mode  centered  at 
0.05pm.  Hence,  the  first  mode  might  well  be  virtually  ignored  in  the  run. 
What  one  needs  to  do  is  to  make  sure  that  quite  a  few  radii  will  be  chosen  to 
span  each  peak,  and  the  simplest  way  to  do  that  is  to  use  the  technique 
suggested  above.  A  more  complicated  way  would  be  to  go  into  the  source  code 
and  change  the  value  of  the  variable  NHALV  to  a  larger  value.  The  number  of 
radii  chosen  between  RLO  and  RHI  if  NBINS  =  0  is  1  +  2**NHALV.  Making  such  a 
change  is  feasible  if  the  kind  of  situation  discussed  is  common  for  a 
particular  user,  but  it  can  exact  a  heavy  time  penalty  by  causing  the  program 
to  utilize  an  unnecessarily  large  number  of  very  large  Mie  size-parameters. 

2.4.3  Subsldaiy  Program  Outputs 

a.  For  all  distribution  types,  the  program  prints  a  quantity 
labelled  "AVERAGE  NUMERICAL  DRY  VOLUME  PER  PARTICLE."  This  quantity  is  used, 
when  appropriate,  to  calculate  the  particle  number  density  DENS. 

b.  For  a  few  distribution  types  the  program  also  prints  a  quantity 
labelled  "AVERAGE  ANALYTIC  DRY  VOLUME  PER  PARTICLE."  This  quantity  can  be 
determined  through  formal  Integration  for  some  of  the  distribution  types  and 
represents  the  result  of  an  analytic  Integration  over  the  radii  between  0  and 
Infinity.  Its  value  may  be  meaningful  in  some  situations:  if,  for  example, 
RLO  and  RHI  were  intended  to  span  virtually  all  particles,  a  comparison  of  the 
"analytic"  and  "numerical"  average  volumes  per  particle  give  at  least  a  useful 
Indicator  of  how  well  that  was  accomplished.  If  the  two  results,  when 
available,  differ  appreciably  it  might  be  advisable  to  think  about  whether  the 
values  of  RLO  and  RHI  are  really  appropriate  for  that  model. 


c.  Hear  the  end  of  the  program,  the  "TOTAL  NUMBER  OF  RADII  USED  HAS" 
is  printed.  The  number  printed  represents  the  sum  of  the  numbers  of  radii  at 
which  Mie  calculations  were  done,  with  all  size-groups  included.  Its  value 
may  be  helpful  in  selecting  a  value  for  DELTA,  but  it  is  printed  mainly  for 
general  interest. 

d.  After  printing  the  values  of  the  phase  functions  (which  is 

omitted  if  parameter  NANG  is  less  than  3),  the  phrase  "TEST  INTEGRAL  OF  PHASE- 
FUNCTION"  followed  by  a  number  is  printed.  If  the  phase  functions  produced  by 
the  run  are  to  be  used  as  input  to  some  subsequent  program  that  will  perform 
numerical  integrations  involving  the  phase  functions,  that  number  should 
always  be  compared  to  unity.  If  its  value  is  appreciably  different  from 

unity,  it  is  inadvisable  to  attempt  to  use  the  phase  functions  in  later 

integrations.  The  discrepancy  is  a  good  measure  of  whether  or  not  enough 
angles  were  used  in  the  run  to  provide  adequate  accuracy  for  numerical 

integrations.  Whenever  such  subsequent  integrations  are  intended,  set  NANG  to 
the  maximum  possible  value  (65  in  the  standard  AGAUS  code).  If  65  angles 
still  do  not  produce  good  agreement  between  the  result  printed  by  the  program 
(that  is,  a  number  very  near  1.0),  a  user  may  be  forced  to  resort  to  use  of  a 
more  sophisticated  numerical  procedure  than  the  trapezoidal  method  used  in 
AGAUS  for  this  test  integration. 

e.  On  occasion,  one  may  find  the  message  "UPPER  MIE  SIZE  LIMIT 

EXCEEDED  IN  INTERVAL  NO."  ...  This  message  indicates  that  at  least  one  value 
of  particle  radius  being  used  for  the  indicated  size-group  was  so  large  that 
the  corresponding  Mie  size-parameter  was  larger  than  ATOP.  This  message  is  in 
the  nature  of  a  warning,  but  there  may  be  no  action  one  can  take  to  alleviate 
the  limitation  unless  ATOP  was  set  below  400.  On  the  other  hand,  if 
sufficient  addressable  memory  is  available,  the  array  A(400)  and  the  parameter 
NDIM  in  subroutine  MIEGX  can  be  increased  substantially.  Experience  with 
MIEGX  has  generally  indicated  that  "numerical  instabilities"  (which  were 
mentioned  in  the  discussion  of  input  data  card  type  3A)  occur  when  the  value 
of  a  size-parameter  exceeds  the  size  of  the  array  A  and  the  parameter  NDIM  of 
MIEGX.  (The  code  attempts  to  compensate  for  such  an  occurrence,  so  the 
program  will  not  ordinarily  abort  as  one  might  suspect.)  Values  of  NDIM  as 
large  as  several  thousand  have  proved  usable  on  systems  that  can  provide  the 
needed  array  space.  The  principal  "problem"  that  has  been  found  with  MIEGX 
for  size-parameters  larger  than  NDIM  has  been  unreliable  values  for  the 
backscatter  efficiency  factor  (and  by  inference,  the  phase  functions  at  large 
angles).  If  those  quantities  are  of  no  interest,  ATOP  can  be  set  to  almost 
any  (large)  value  and  the  extinction  and  scattering  cross  sections  may  be 
quite  reasonable.  Users  can  test  this  for  themselves  by  exercising  the  MQRTE 
=  12345  option  (or  by  running  model  7)  and  examining  the  behavior  of  the 
efficiency  factors  as  a  function  of  radius  or  size-parameter.  The  extinction 
and  scattering  efficiency  factors  should  behave  smoothly  and  approach  the 
value  of  2.0  for  very  large  size-parameters.  The  smooth  behavior  is  also 
expected  for  the  backscatter  (radar)  efficiency  factors  for  large  size- 
parameters.  A  look  at  the  values  of  the  latter  may  demonstrate  the  "problem" 
under  discussion  here. 

f.  Certain  anticipated  error  conditions  also  generate  output 
statements  aimed  at  making  the  source  of  error  reasonably  visible  to  the 
user.  In  most  cases  these  errors  are  deemed  "fatal"  and  the  run  is 
terminated. 


2.5  The  WAVE  Looping  Options 


AGAUS  permits  looping  over  any  parameter  that  appears  on  data  card  type  5  (see 
section  2.3.5)  if  an  appropriate  input  data  set  is  prepared.  This  section 
discusses  one  possible  application  of  this  option.  Suppose  that  one  wishes  to 
examine  the  way  in  which  the  extinction  coefficient  of  a  particular 
hygroscopic  aerosol  model  changes  with  relative  humidity,  and  that  data  for 
EMUA  (17(f))  are  available  for  five  values  of  the  relative  humidity.  One  would 
then  set  parameter  NWAVE  =  5  on  data  card  number  2  and  prepare  cards  of  types 
3  and  4  in  the  usual  way.  The  choice  NWAVE  =  5  means  that  at  least  five  of 
the  type  5  data  cards  are  needed--one  for  each  value  of  the  relative  humidity 
(RELHUM).  These  five  cards  would  ordinarily  carry  identical  data  for 
parameters  WAVE,  EMA,  CAYA,  RHOA,  CONC,  and  TEMP,  but  different  values  of 
RELHUM  and  EMUA.  The  value  of  EMUA  on  a  given  card  should,  of  course, 
correspond  to  the  value  of  RELHUM  on  that  card. 

The  use  of  five  type  5  cards  to  handle  five  values  of  RELHUM  assumes  that  only 
one  size-distribution  is  in  use  (that  is,  that  NIDSTP  =1).  If  NIDSTP  is 
greater  than  1,  the  required  number  of  type  5  cards  is  given  by  the  product  of 
NWAVE  and  NIDSTP.  As  an  example,  suppose  that  NWAVE  -  5  and  NIDSTP  =  2. 

There  must  then  be  ten  of  the  type  5  data  cards  arranged  as  follows: 

WAVE,  EMAj,  CAYAlt  RHOAj ,  CONCp  RELHUM^  TEMP,  EMUAn 

WAVE,  EMA2,  CAYA2,  RH0A2,  C0NC2,  RELHUMj  TEMP,  EMUA12 

WAVE,  EMAj,  CAYAlt  RHOAj,  C0NClf  RELHUM2,  TEMP,  EMUA21 

wave,  ema2,  caya2,  rhoa2,  conc2,  relhum2,  temp,  emua22 


WAVE,  EMAl,  CAYAlt  RH0Alf  CONCj,  RELHUMj  TEMP,  EMUAgj 
WAVE,  EMA2,  CAYA2,  RH0A2,  C0NC2,  RELHUMg,  TEMP,  EMUA52 

where 

EMA^ ,  CAYAj  are  the  optical  constants  for  (dry)  aerosol  component 
1  at  wavc’ength  WAVE  and  temperature  TEMP, 

RH0A.j ,  CONC^  are  the  mass  density  and  mass  concentration  for 
component  i , 

RELHUMj  is  the  jth  value  of  relative  humidity, 

and 

EMUA44  is  the  value  of  TT(f )  for  aerosol  component  i  at 

relative  humidity  value  RELHUMj. 


In  this  example,  parameters  WAVE,  EMA.,* ,  CAYAi ,  RHOA^ ,  CONCi  and  TEMP  are  held 

constant,  with  only  RElHUMj  and  EMUA,-,;  changing  from  card  to  card.  One  could 

set  up  a  data  deck  for  looping  over  any  other  parameter  on  the  type  5  cards  in 
a  similar  fashion  if  desired:  use  (NWAVE)  x  (NIDSTP)  type  5  cards  with  all 
parameters  fixed  except  the  one  over  which  looping  is  to  occur.  Note, 
however,  that  AGAUS  does  not  "remember”  any  of  the  parameters,  so  valid  values 
of  WAVE,  EMA,  etc.,  mist  appear  on  every  individual  card.  There  is  only  one 
exception  to  this  rule:  If  IW  =  0  (water  only  case),  then  the  values  of  EMA 
and  CAYA  will  be  found  from  the  internal  tables;  any  input  value,  including 
blanks,  will  be  ignored. 

2.6  Differences  between  AGAUS  82  and  Previous  Versions 

The  structure  of  AGAUS  82  differs  in  several  ways  from  that  of  the  version 
listed  in  reference  9.  The  principal  differences  are  summarized  here  for  the 
possible  benefit  of  users  of  the  previous  version. 

a.  The  ability  to  "loop"  over  a  mixture  of  size-distribution  models 

has  replaced  the  loop  over  several  components  (of  varying  optical  properties) 

having  a  fixed  size-distribution  model.  Each  component  now  requires 

definition  of  a  separate  size-distribution  model  even  if  several  components 

are  to  be  given  the  same  size-distribution  model.  This  option  is  not 
available  for  the  "arbitrary"  (IDSTP  =  0)  model. 

b.  Several  new  output  quantities  have  been  added  with  breakdowns 

into  particle  size  groups.  These  Include  quantities  such  as  the  number  of 
particles  per  gram  and  cross  sections  in  square  centimeters  per  gram. 

c.  The  normalization  of  the  phase  functions  has  been  changed  so  that 
their  integral  over  all  solid  angle  should  yield  unity. 

d.  The  algorithms  for  handling  the  growth  and  variations  in  optical 
properties  of  hygroscopic  aerosols  as  functions  of  relative  humidity  have  been 
replaced  by  somewhat  more  sophisticated  ones. 

e.  The  original  Mle  subroutine  MIEGX  has  been  replaced  by  one  using 
different  numerical  methods  and  is  valid  for  larger  size-parameters  than  the 
previous  version. 

f.  Unweighted  integrations  of  extinction  coefficients,  phase 
functions,  etc.,  over  wavelength  have  been  deleted,  and  the  limit  of  10 
wavelengths  per  run  has  also  been  deleted. 

g.  Some  of  the  IEO  options  have  been  omitted. 

h.  All  Intrinsic  "bimodal"  size-distribution  models  have  been 
deleted,  as  has  the  one  labelled  "Del rmendjian's  Model.”  Bimodal  cases  can 
now  be  treated  by  using  the  looping  option  mentioned  as  item  a  above.  A  new 


9D.  Delrmendjian,  1969,  Electromagnetic  Scattering  on  Spherical 
Polydispersions,  American  Elsevier  Publishing  Company,  Inc.,  New  York 
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nybrid  lognormal /power- law  size-distribution  model  das  been  added,  and  a 
special  "model"  ( I DST?  =  7)  has  been  added  to  print  tabulations  of  Mie 
efficiency  factors. 

i.  Tne  option  for  generating  Gauss-Legendre  expansion  coefficients 
for  the  phase  function  has  been  removed  as  has  the  use  of  special  angles 
associated  with  that  option. 

j.  A  new  input  paraneter  (ATOP)  allows  users  to  set  an  upper  limit 
on  the  .lie  si ze-parameters  to  De  used  in  a  run. 

k.  The  ordering  of  the  input  data  cards  and  the  data  belonging  to 
eacn  "type"  of  card  have  been  changed  substantially. 

l.  A  new  subroutine  (DISTR)  handles  the  calculation  of  the 
distribution  functions  f ( r)  for  all  cases  except  the  arbitrary  (IOSTP  =  0) 
case.  Use  of  tnis  routine  makes  possible  removal  of  the  upper  limit  of  513 
radii  per  wavelength  (except  for  IUSTP  =  0  cases). 
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